% x - x values where predictions are made
% m - mean, should be same length as x
% sd - standard deviation, sould be same length as x
% obs_x - x values of observations
% obs_y - y values of observations, should be same length as obs_x
% bounds - [min_x max_x min_y max_y] bounds for axes
% dot_size - size of observation dots
% x_label - x label text
% y_label - y label text
% obs_label - text to use for observation dots in legend
% legend_location - where to place legend (e.g. 'NorthEast')
% width - width in centimeters
% height - height in centimeters
% name - filename to save plot to
% legend_location2 (optional) - if you'd prefer to have two
%   legends, one for observations and the other for mean/sd, the
%   location of the mean/sd legend.

function make_gp_plot (x, m, sd, obs_x, obs_y, bounds, dot_size, ...
											 x_label, y_label, obs_label, legend_location, width, ...
											 height, name, legend_location2)

split_legend = (nargin > 14);

close all

x = x(:);
m = m(:);
sd = sd(:);
sd(isnan(sd)) = 0;
sd(sd == 0) = eps;
m(isnan(m)) = 0;

hold on;
SDh = fill([x; x(end:-1:1)], ...
  [m + 2 * sd; m(end:-1:1) - 2 * sd(end:-1:1)], ...
  [0.87 0.89 1], 'EdgeColor', 'none');
meanh = plot(x, m);
observationsh = plot(obs_x, obs_y, '.');

fh = figure(1);

set(fh, 'units', 'centimeters', ... 
  'NumberTitle', 'off', 'Name', 'plot');
pos = get(fh, 'position'); 
set(fh, 'position', [pos(1:2), width, height]); 

axis(bounds)

set(gca, 'TickDir', 'out')
set(gca, 'Box', 'off', 'FontSize', 10); 
set(fh, 'color', 'white'); 
set(gca, 'YGrid', 'off');

set(meanh, ... 
  'LineStyle', '-', ...
  'LineWidth', 0.75, ...
  'Marker', 'none', ...
  'MarkerSize', 10, ...
  'Color', [0 0 0.8] ...
  );

set(observationsh, ... 
  'LineStyle', 'none', ...
  'LineWidth', 0.5, ...
  'Marker', '.', ...
  'MarkerSize', dot_size, ...
  'Color', [0.2 0.2 0.2] ...
  );

xlabel(x_label)
ylabel(y_label)

ah = gca;

if (split_legend)
    hLegend = legend(...
        ah, ...
        observationsh, ...
        obs_label, ...
        'Location', ...
        legend_location, ...
        'boxoff' ...
        );
    
    legend boxoff
    
    ah2 = axes('position', ...
               get(gca, 'position'), ...
               'visible','off');
    
    hLegend2 = legend(...
        ah2, ...
        [meanh, SDh], ...
        'mean' , ...
        '$\pm 2$ SD$', ...
        'Location', ...
        legend_location2 ...
        );
    
    legend boxoff
    axes(ah)
else
    hLegend = legend( ...
        [observationsh, meanh, SDh], ...
        obs_label, ...
        'mean' , ...
        '$\pm 2$ SD', ...
        'Location', ...
        legend_location, ...
        'boxoff' ...
        );
    
    legend boxoff
end

set(0, 'defaulttextinterpreter', 'none')
matlabfrag(name);
